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A new Direct Numerical Simulation (DNS) of 
Couette-Poiseuille flow at a higher Reynolds number 
is presented and compared with DNS of other wall- 
bounded flows. It is analyzed in terms of testing semi- 
theoretical proposals for universal behavior of the ve- 
locity, mixing length, or eddy viscosity in pressure 
gradients, and in terms of assessing the accuracy of 
two turbulence models. These models are used in 
two modes, the traditional one with only a dependence 
on the wall-normal coordinate y, and a newer one in 
which a lateral dependence on z is added. For pure 
Couette flow and the Couette-Poiseuille case consid- 
ered here, this ^-dependence allows some models to 
generate steady streamwise vortices, which generally 
improves the agreement with DNS and experiment. 
On the other hand, it complicates the comparison be- 
tween DNS and models. 


1 Introduction 

The present study extends that of Johnstone, Cole- 
man and Spalart in 2010, in which a single case of 
Couette-Poiseuille flow at a relatively low Reynolds 
number was presented. The purpose was to extend 
DNS conducted between two parallel plane walls to 
flows with a range of pressure gradients, relative to the 
common case of Poiseuille flow which has only a weak 
favorable gradient, without facing the difficulties of a 
spatially developing boundary layer. Compared with 
Poiseuille flow, Couette-Poiseuille requires very little 
extra coding, but raises severe issues of domain size 
because of the global vortices. The earlier study ap- 
peared to provide a fairly clear if surprising answer 
to the seemingly simple question of how the law or 
laws of the wall extendfs) to flows in pressure gradi- 
ents. The question was couched as follows. 

Assuming a constant-stress ‘Couette-flow’ region 
in wall-bounded turbulence leads to at least three clas- 
sical relationships between the derivative of the mean 
velocity dU /Ay and the shear stress —u'v', the two 
quantities at the core of turbulence modeling. Ex- 
pressed in terms of length scales to allow direct com- 


parisons, these include the one leading to the logarith- 
mic law of the wall for the velocity. 


u T 

AU/dy = KV ' 


( 1 ) 


where n is the Karman constant, u T = \Jt w / p 
the surface-friction velocity and y w the wall-normal 
distance. In terms of eddy viscosity v t = 
— u'v' / (dU/dy), we also have 


v t —u'v'/u T 
^ = dU/dy = KVw 


( 2 ) 


while consideration of the mixing length £ = 

(v t /(dU/Ay)) 1/2 gives 


V 

= -WJaH = Ky ’ 


—u'v' 


(3) 


Since —u'v' = r(y w )/p = u/. in a constant-stress 
layer and outside the viscous region, all three expres- 
sions are equivalent, and arguing between them proba- 
bly has limited value. None of the three has a system- 
atic derivation. With pressure gradient, r/p deviates 
from ui and the three proposals conflict, so that at best 
one of the three will remain effective. This could make 
such a flow very discriminating, while very accessible 
to DNS. Ideally, in the near future DNS will be con- 
ducted over a fair range of Couette-Poiseuille cases, 
and new experiments with very fine measurement of 
the skin friction will be conducted for support and to 
reach more convincing Reynolds numbers. 

Our early results (Johnstone et al. 2010) indicated 
that the logarithmic -law length scale (1) is much less 
affected by pressure gradients than the other two op- 
tions, a finding for which we had (and continue to 
have) no theoretical justification, nor any strong sense 
of whether it might be exact, or is only approximate. 
The demonstration had a preliminary character, due 
to Reynolds-number limitations and there being only 
one case. The thickness of the inertial sublayer es- 
pecially on the adverse-gradient (APG) side was not 



large enough to convince all observers. We address, 
to an extent, this shortcoming in the current computa- 
tions. 

This study can also be viewed as a computational 
counterpart to the boundary-layer experiments per- 
formed by Galbraith & Head (1975) and Galbraith, 
Sjolander & Head (1977), who long ago were also 
interested in the relative sensitivity to pressure gradi- 
ents of the law of the wall (1), eddy viscosity (2) and 
mixing length (3). Our conclusions again will turn 
out to be broadly compatible with theirs. We expect 
boundary-layer DNS to produce the same conclusion, 
which was already hinted at in Spalart & Watmuff’s 
(1993) adverse-gradient boundary-layer DNS. 

Another new element of the present study is results 
of Reynolds-Averaged Navier-Stokes (RANS) mod- 
els. These are exercised with and without streamwise 
vortices, a mode of operation recently proposed with 
some success by Spalart, Garbaruk & Strelets (2014) 
for pure Couette flow. Such vortices need to be con- 
sidered if a rigorous comparison of DNS and RANS 
results is to be performed, and it is possible their pres- 
ence also alters the answer to the seemingly simple 
question regarding the alternative length scales formu- 
lated above. 

2 Approach and Problem Definition 

In order to introduce favorable and adverse 
pressure gradients and to separate them from 
the Reynolds number, we consider incompressible 
Couette-Poiseuille flow, driven by a steady uniform 
pressure gradient d(P/p)/dx and in-plane stream- 
wise wall motion, in opposite directions, of speed 
U w . This is a two-parameter family of flows. A 
Reynolds number, based on the wall-velocity differ- 
ence AU = ‘211 w and full channel height 2 h, of 
Re = 80, 000 is prescribed (i.e. twice that of our 
2010 study), along with a streamwise pressure gradi- 
ent d(P/p)/dx = —0.0014(7j/ft.. This leads to non- 
dimensional pressure gradients p + = dr + /dy + = 
[d(P/p)/dx\ [v/ul\ of p+ « -0.0003 and +0.0030 
on the favorable-pressure-gradient (FPG) and adverse- 
pressure-gradient (APG) sides, respectively. The 
FPG value of p + is the same as in Poiseuille flow 
with Re r = u T h/v = 3125, and therefore quite 
weak. The simulation uses a parallel version of 
the Fourier/Chebyshev spectral channel code of Kim, 
Moin & Moser (1987), from which it differs algorith- 
mically only in that a third-order Runge-Kutta/Crank- 
Nicolson scheme is used for the time integration 
(which is formulated in the reference frame moving 
at the average speed of the two walls). 

The simulation was initialized by random distur- 
bances (of amplitude large enough to trigger non- 
linear/bypass transition) on an under-resolved grid at 
lower Reynolds number, and continued until the flow 
reached a statistically stationary state, after which the 
Reynolds number was increased to its final value and 


Table 1: Run parameters. 


Wall 

U T /^Uuj 

p+ 

Ax + 

Az+ 

Vio 

APG 

0.0145 

+0.00285 

5.4 

2.7 

3.0 

FPG 

0.0302 

-0.00032 

11.3 

5.6 

6.1 


the flow re-projected onto finer grids multiple times, 
each advanced to a stationary condition. The results 
shown below are from the fully resolved state, run suf- 
ficiently long (see below) to obtain well-defined statis- 
tics averaged in time and in the two homogeneous di- 
rections. 

At this Reynolds number, a real-space resolution 
of n x x n y x n z = 1,344 x 281 x 1,344 colloca- 
tion/quadrature points is required to capture all rele- 
vant scales throughout the domain. (The 3/2 rule is 
applied in x and z for dealiasing.) This grid count, 
near 0.5 x 10 9 , does not allow multiple simulations to 
be performed easily. In wall units from the FPG side 
(i.e. the side with larger u T ), this corresponds to a time 
step A/ + ss 0.24, streamwise and spanwise grid spac- 
ing of A.r + ss 11 and Az + ~ 6, and wall-normal 
distance of the 10th grid point y+ 0 ss 6. This is more 
stringent than the ‘conventional’ guidelines for a spec- 
tral method, of (Ax + , y+ 0 , Az + ) < (20, 10, 7), which 
is preferable when seeking definitive results from DNS 
(Spalart et al. 2009). 

The size of the computational domain is 47 rh x 
2 h x 27 ih, in the streamwise x, wall-normal y and 
spanwise z directions, respectively. The spanwise do- 
main length L z is large enough to contain two pairs 
of large-scale streamwise roll cells, which tend to 
reach from wall to wall (Figures 2 and 3), typical of 
those found in pure Couette flow. Johnstone (2014, 
personal communication) observed similar cells in 
his lower-Re, larger-domain Couette-Poiseuille DNS, 
for which dU w h/v = 40,000 and d{P/p)/dx = 
— 0.00085(7^,//i, with (L X1 L Z ) = (167t/i, dirh). It is 
very likely that much stronger pressure gradients, pro- 
ducing flows closer to the Poiseuille limit than those 
considered here, are required to prevent the formation 
of large-scale, domain-filling streamwise rolls cells. 
The presence of the rolls complicates the interpreta- 
tion of this flow as a prototype of canonical FPG and 
APG boundary layers, and thus care must be exercised 
when using the present results to draw general con- 
clusions regarding turbulence theory and modeling. In 
particular, the results may depend on the size of the 
spanwise period, which is arbitrary. This issue will be 
revisited below. 

3 Results 

Overall character of the flow 

Figure 1 illustrates the asymmetry in the velocity 
profile introduced by the pressure gradient and mov- 
ing walls, which still has a strong Couette character 
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Figure 1: Top, mean velocity; bottom, Reynolds and total 
shear stress. 
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Figure 2: Contours of instantaneous streamwise velocity in 
vertical-spanwise ( y-z ) plane from realization at 
(top) beginning and (bottom) end of averaging pe- 
riod. 


ii/Uw: -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 



z/h 


Figure 3: Contours of streamwise-averaged streamwise ve- 
locity in vertical-spanwise (y-z) plane from real- 
ization at (top) beginning and (bottom) end of av- 
eraging period. 


with only one sign of dU / Ay and r. The effective FPG 
and APG wall layers are characterized by the negative 
and positive (with respect to the reference frame at- 
tached to each moving wall) wall-normal total shear- 
stress gradients dr/dy w also shown in Figure 1. The 
level of equilibrium of the time- and ( x , z)-plane av- 
erage, which was gathered over a period of approxi- 
mately 460/i/C/u,, can also be inferred from this figure 
by comparing the total stress profile to the linear ide- 
alization defined by the pressure gradient. This figure 
also shows that the viscous shear stress is confined to 
thin regions, compared with h; this widens the region 
over which pure turbulent behavior can be hoped for. 

Although the linear stress profile reflects a mature 
solution very close to equilibrium, it appears the roll- 
cell structures may not have reached their final (sta- 
tionary) configuration at the beginning of the averag- 
ing period. This is suggested by the difference be- 
tween the streamwise-velocity contours in y-z planes 
at times corresponding to the beginning (top figures) 
and end (bottom figures) of the averaging, which are 
shown in figures 2 and 3 (the former showing a sin- 
gle y-z plane, the latter the streamwise-domain aver- 
age). The difference is consistent with the roll struc- 
ture ‘splitting’ laterally from one to two, transitioning 
from a state in which the roll energy resides primar- 
ily in the (k* x ,k* z ) = (k x L x /2Tr 1 k z L z /2n) = (0,±1) 
mode (where k x and k z are the streamwise and span- 
wise wavenumbers) to one in which (0, ±2) domi- 
nates. This has been confirmed by examination of the 
velocity spectra from the channel centerline at the two 
times in question. The DNS will need to be contin- 
ued to ascertain whether the roll’s final state remains 
‘locked’ into the (0, ±2) mode, and to quantify the ef- 
fect the roll evolution has on the mean profiles shown 
below. 
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Figure 4: Top, Reynolds- and total-shear stress; bottom, 
root-mean-squared (RMS) velocity fluctuations, in 
local wall units. 


Figure 4 initiates the comparison of the two wall 
regions, each with a normalisation based on its own 
friction velocity. The influence of the pressure gradi- 
ent on the near-wall variation of the Reynolds stresses 
is very visible. The thick-solid line near y + = 0 gives 
the slope, in wall units, of the total shear stress. The 
pressure-gradient effect reaches all the way to the wall, 
especially for the two wall-parallel stresses u'u' and 
w'w'\ in other words, they grossly violate any Taw of 
the wall,’ defined solely in terms of y + and u T . Note 
however that even without pressure gradient, when 
varying Reynolds number the diagonal stresses violate 
what would be a law of the wall to a similar degree, 
at least for Reynolds numbers currently accessible to 
DNS (Spalart 1988, Hoyas & Jimenez 2006). 



Figure 5: Turbulence length scales. 


Comparison between length scales 

In Figure 5 we re-visit the question of the pressure- 
gradient response of the wall-turbulence semi-theories 
(1), (2) and (3). As in Johnstone et al. (2010), the log- 
law version of the length scale (1) is found to be sig- 
nificantly less vulnerable to pressure gradient than the 
other two. This is seen more clearly with the higher 
Reynolds number, and applies to both walls. Note that 
the line l = ny w is not in-between (1) and another 
curve; in other words, the results do not indicate a 
mere preference towards (1), they suggest that (3) and 
especially (2) grossly fail to match reality. This argues 
against the success of a theory based on a weighted 
average of the three candidates. 

Compared with our first Couette-Poiseuille paper 
(Johnstone et al. 2010), the higher Reynolds number 
allows us to refine the conclusions: especially in the 
region below roughly y/h = —0.7, length scale (1) ap- 
pears to have a smooth and consistent deviation from 
KoUw, where kq = 0.4. One interpretation could be 
that this region has a lower value of k, of the order 
of 0.375; this value is quite low, even with the current 


uncertainty over the exact value of n. Also note that 
the same attitude applied near y = +0.8 and therefore 
on the APG side would suggest a higher value than 
0.4 for k. Another interpretation is that the Reynolds 
number is not sufficient for the result to arrive at the fi- 
nal behavior. However, considering that the Reynolds 
number has been doubled, it now appears established 
that the ‘resilience’ of (1) we noted in our earlier study 
is not perfect. 

The trend seen here is consistent with Galbraith’s 
1975 experimental findings, which we believe have not 
been adequately recognized, confirmed, challenged, or 
adopted as a guide for turbulence modeling. Unfortu- 
nately, this line of work remains narrow and isolated. 

Detailed behavior of the velocity profile 

Having established that the velocity profile is 
closer to universal behavior than the other candidates, 
we now study it more closely. The effect of the pres- 
sure gradient on the law-of-the-wall scaling of the 
mean velocity is shown in Figure 6, which includes 
results from both the present and earlier Couette- 
Poiseuille DNS, along with the pure-Poiseuille pro- 
file from Hoyas & Jimenez (2006), for which p + = 
—5.0 x 10 -4 (indicating a weak FPG). The higher 
Reynolds number of the new simulation is evident in 
the clearer distinction between the buffer, (nominally) 
logarithmic, and wake regions for both the FPG and 
especially the APG sides. The FPG profile agrees well 
with the reference higher-Reynolds -number profile of 
Hoyas & Jimenez, while the APG result falls well be- 
low the channel benchmark: the viscous sublayers all 
agree to plotting accuracy, but at larger y + the APG 
profiles dip below the FPG mean; the deviation be- 
comes noticeable near y + = 15. Note that the grid 
resolution in wall units is more favorable for this wall, 
with its lower u T , so that a higher level of numerical 
error near this wall is not suspected. At y + = 30, the 
total shear stress is about 8% higher than at the wall, 
which can already be considered a significant devia- 
tion. In other words, if measured by the shear stress the 
FPG cases are actually clustered very close together, 
and the APG case is isolated; this makes the behavior 
of the velocity less surprising. 

A more sensitive measure of the velocity profile is 
presented in Figure 7, in terms of what we refer to as 
the Karman measure km = d(lnr/ + )/dt/ + . It is equal 
to the length scale in (1), divided by y + . Also included 
are the Hoyas & Jimenez (2006) channel data and re- 
sults from a series of Ekman-layer DNS over a range 
of Reynolds numbers (Spalart et al. 2009; Johnstone 
2012). The Karman measure equals the Karman con- 
stant k in an exactly logarithmic region; the lack of a 
plateau (expected to begin around y + = 100, roughly) 
indicates that such a layer is not present for the flows 
examined here. This may be blamed on the relatively 
low Reynolds numbers amenable to DNS, but perhaps 
not totally, considering that values more than 20 times 
larger than in the pioneering DNS of the 1980s have 



Figure 6: Velocity profiles in wall units. The plane 

Poiseuille results of Hoyas & Jimenez (2006) are 
included for comparison (p + = —5.0 x 10~ 4 ). 



Figure 7: Karman measure. Included for comparison are 
DNS results from the plane-Poiseuille flow of 
Hoyas & Jimenez (2006) and the Ekman-layer 
profiles of Spalart et al. (2009) and Johnstone 
( 2012 ). 


been reached. If the current thinking among some ex- 
perimentalists (based on very high Reynolds number 
pipe-flow measurements) - that a constant-K region 
should only be expected for Reynolds numbers large 
enough that y + = 10 3 falls well below 0. 1(5 (where 5 
is the layer thickness) - then a ‘DNS answer’ for the 
value of k - and confirmation that it is indeed universal 
- is many years away. 

For now, we note that Figures 6 and 7 support 
the validity of the general law of the wall (i.e. that 
dU + /dy + depends only on y + ), which is unaffected 
by system rotation (acting in the Ekman flow) and the 
sign and magnitude of pressure gradient (in the pure- 
and Couette-Poiseuille channels), for 0 < y+ < 35. 
Focusing only on the FPG results, a further case can 
be made that the range of established universality ex- 
tends to y + ss 400 or so, since the rotating-boundary- 
layer and channel results agree here. The trend versus 
Reynolds number is very encouraging. Although some 
questions remain open, this finding tends to contradict 
claims made by some that even the high-Reynolds- 
number limit of k would be flow-type dependent, that 
is, different for example in channels, pipes and various 
versions of boundary layers. The behavior in APG is 
a little dissonant, but we pointed out that its pressure 
gradient in wall units is still fairly strong. 
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Figure 8: Values of U + at y + = 50 for various pressure 
gradients, from experiment and simulation. New 
results are shown by solid red upward triangles. 


We return to the pronounced downward shift with 
increasing p + seen for y + > 15 in Figure 6. Fig- 
ure 8 (essentially the same as Figure 6 in Johnstone 
et al.) documents the variation of (7 + (50), the ve- 
locity at y + = 50, for a range of pressure gradients. 
This includes a wide group of earlier DNS and experi- 
ments, some with much stronger values of p + than the 
present ones, in both directions. The new results are 
in line with the old ones. Both our old and new APG 
results are somewhat below the general trend. The re- 


sults are compatible with the conjecture that U + ( 50) 
primarily depends on p + , as opposed to other param- 
eters such as Reynolds number. The general conjec- 
ture would be that the velocity profiles follow a law 
of the type U + = f(y + ,p + ), with the flow Reynolds 
number and other measures of the outer flow having a 
much weaker effect than p + does, over a useful range 
of y/6 where 5 is a measure of the outer-layer thick- 
ness. Such a finding could be used to construct a 
pressure-gradient-dependent wall function for U + in 
non-separated boundary layers. However, there are 
many issues with this, such as the difference between 
p + and dr + /dy + in flows other than channels (John- 
stone et al. 2010). Finally, recall that mixing-length 
reasoning would predict a positive slope in this fig- 
ure (Spalart & Watmuff 1993), thus repeating the egre- 
gious failure of that approximation seen in Figure 5. 




Figure 9: Karman measure implied by SST and SA models. 


Results of turbulence models 

We conclude by examining the ability of two 
widely-used RANS models to reproduce the DNS re- 
sults, with attention given to possible roll-cell effects. 
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Figure 10: Turbulence length scales implied by SST and SA 
models. 


The models are Menter’s (1994) SST and the Spalart- 
Allmaras (1994) (SA) models, the latter with and with- 
out the quadratic-constitutive-relation (QCR) option 
(Spalart 2000). The SA/QCR model is exercised in 
the spanwise-periodic mode introduced, in a study of 
pure Couette flow by Spalart et al. (2014), with vari- 
ous values for the period L z , which controls the ‘roll 
cells.’ Briefly, the approach is to obtain RANS solu- 
tions that are steady and independent of x, but have a z 
dependence and three velocity components: ((7, V, W) 
are functions of (y, z). Spalart et al. (2014) found that 
eddy-viscosity models combined with a simple QCR 
(and only these models, so far) sustain roll cells that re- 
semble those known in Couette-flow experiments and 
DNS (evoking Figure 3), and account for sometimes 
half of the momentum transfer in the core region of 
the flow. It was also found that this transfer depends 
quite strongly on the rolls’ spanwise spacing, which is 
arbitrary within a range; that is, the equations accept 
the rolls only over a band of wavelengths. Note that 
the ( k *, k*) = (0, ±2) mode associated with the DNS 
rolls in Figure 3(bottom) implies their lateral spacing 
is A- = L z /2 = nh, such that an SA/QCR run with 
L z /h = 7 r corresponds to the most mature DNS spac- 
ing. 

The Karman-measure profiles are shown in Fig- 
ure 9. The SST model is not very accurate in the buffer 
layer; this is a known issue, ultimately tied to the non- 
use of a viscous damping function similar to f v i in the 
SA model. Farther from the wall, SST is in noticeably 
better agreement with DNS than SA, unless SA/QCR 
is used with the most favorable value of L z /h, namely 
2 (a value which gives the best results for both walls, 
but is not very close to the current DNS value, namely 
7r). This confirms that the roll cells have a major in- 
fluence on the RANS results, not only in Couette flow 
but also in the present blend with Poiseuille flow. Runs 
with SST-QCR and roll cells will also be instructive, 
although the use of the model without roll cells is cer- 
tainly appropriate. 

In Figure 10, we examine the three length scales 
as produced by the turbulence models. It confirms 
Figure 9 in terms of ranking the models for accuracy. 
The figure also gives a commentary on the nature of 
transport-equation turbulence models, as opposed to 
algebraic models. A mixing-length algebraic model 
would agree precisely with equation (3), for instance, 
up to some distance from the wall of order h. This does 
not happen with either of the present models. These 
models were constrained to agree with (1), (2) and (3) 
in a constant-stress layer, and it could have been that 
they fortuitously agree with one of the three in general. 
Instead, they agree with (1) better than (2) or (3); this 
is a favorable trend, since this is also what DNS does, 
but this is not by design and simple behavior is not to 
be expected. In the SA model, the behavior could be- 
come ‘by design’ if the f w function were adjusted to 
match new results. 
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4 Summary and Future Work 

A new Couette-Poiseuille DNS is used to explore 
the universality and robustness of conflicting versions 
of near-wall similarity. The impact of pressure gra- 
dient on the log-law relationship for the mean veloc- 
ity is found to be significantly weaker than it is on 
eddy-viscosity and mixing-length alternatives. Com- 
pared to our previous findings, the present conclusions 
are less clouded by low-Reynolds-number effects, and 
are thus more compelling. On the other hand, the ef- 
fect of the large-scale streamwise roll cells/vortices 
suspected in the DNS is not yet fully resolved, and 
we must not over-estimate the authority of RANS re- 
sults for such phenomena. The reality is that Couette- 
Poiseuille DNS is far from immune to the problem of 
domain size that affects Couette flow. 

We again take note of an encouraging trend among 
many datasets for a possible pressure-gradient-based 
correction of the log law in Figure 8, but recognize 
that a much larger body of evidence is needed before 
serious proposals are made. 

Further DNS (now underway) is needed to com- 
plete the study. We aim to quantify the evolution of 
the roll structure, and its effect on both the dynamics 
of the flow (in terms of momentum transport) and the 
implications for turbulence theory and RANS mod- 
eling. It is remotely possible that the ranking of the 
‘winners’ in the length-scale ‘contest’ (Figure 5) will 
be re-ordered once this is done. This will involve ex- 
amining time- as well as streamwise-averaged fields, 
to ascertain whether a cell remains at a fixed spanwise 
location, meanders in time and/or streamwise direc- 
tion, or even has a finite life. A plausible separation of 
the Reynolds stresses into ‘roll-cell’ and ‘small-scale’ 
contributions would be used after studying the two- 
dimensional spectra in the x-z plane. Statistics gath- 
ered from fields that have had the energy of the modes 
corresponding to the cells removed would be used to 
recompute the length-scale profiles examined above. 

In the longer term, the pressure gradient and the 
wall-velocity difference will not have to be in the same 
direction, so that skewed wall-bounded flows can be 
simulated by the same method, and compared with the 
Ekman layer. 
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